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1. The Interpolation Problem 


Given the set of data points (xp ¥p oF) » it is desired to construct a 
function F(x,y) , such that F(X) s¥ 4) = Fok = 1,...,n . For large sets 
of data it is desirable for the method to be local, that is, the value of 
F(x,y) depends only on the value of f, at nearby points (xy 5¥p) ; 

This problem is receiving a great deal of attention and discussions of 
it and proposed methods can be found in [1], [2], [4]. [7]. and [8]. 

2.0. The Interpolation Scheme 
This is a local method, the general idea having been discussed in [4]. 


The basic idea is to construct local interpolants, F. , which are then 


g 
weighted by functions W, having limited support to obtain the function 


(1) F(xsy) = 2 Wi (xsy)Fo(xsy)/) Wy (x.y) 
L & 


The details are fully discussed in the reference, but the important fact is 
that F(x,y) will take on the value f. at (x) s¥,) if for each 2% where 
W (xp s¥y) =z0, Fo (xp yy) = f, . In the referenced paper, the weight func- 
tions Wo were taken to be of the form 


- 3(d,/r,)° + 2(d,/r,)° ee 


QL L 
Wo (xsy) = 
0 » d>r, 
where vr, is the radius of the smallest circle centered at (x) »¥p) which 


contains a given fixed number of data points, and d, is the distance from 
(x,y) to (x, ¥y) ; 

This scheme, used with local interpolants, F. which were taken to be 
certain optimal approximations, yielded reasonably good results. However, the 
computational burden was rather high. This was due in part to a great deal of 


overlap in the regions where the W, are nonzero. In addition, the approxi- 
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mation is defined only on the union of the circles around each of the data 
points, which could cause problems. 

One of the advantages of using optimal approximations is that the 
basis functions are generated by the (x, 5¥p) and the system is automat- 
ically nonsingular. Within some limitations, the resulting equations can 
easily be solved by Cholesky decomposition [5]. This opens the way for the 
present approach, which is to choose rectangular regions on which weight 
functions are non-zero, thus being able to carefully govern the amount of 
overlap with a resulting decrease in the necessary computations, as well as 
simplification of the weight functions. 

2.1. Weight Functions 

With these ideas in mind, we are now ready to describe the present 

selection of regions over which the weight functions are non-zero. These 


regions will be rectangles defined by the following parameters. Let ny 


~ 


and n, be given positive integers and let finite values of Ky oXo aves aXp 
x 


and Yq Vor-++ Vp be given. For convenience in this section we let 


~ ~ 


ae a and x 


‘ ¥ ny =o. For each i= 1,....n, and 


j= Taeeeomy » let Ris denote the rectangle (X54 X44) X51 jar) , 


+ 
ny ] 


Let He(s) =]- s3(6s" - 15s + 10) , the Hermite quintic satisfying 
H,(0) =1, Hg" (0) = He" (0) = He (1) = He' (1) = He" (1) = 0. We then define 
functions which are piecewise quintics with continuous second derivatives 
which have the property that they are non-zero on two intervals and satisfy 
» j= 0,1,....n, +1 


eee > i= stgeeced + ] . 
on i= 0,1 ny 


In particular, 


] me PY X< x 
xX ad Xy 
V(x) = H5(e = ) 2X) Sx < % 
2 |. 
0 .s Xo 
0 aX < Xe y 
V(x) = x-X 
j i ~ 
He ( = 2 ) Py Xj sX< Xia] 
i+] j 
0 >» X2 Xoa] 
for i = 2,.. a, 1, and 
0 +x <xX 
a 
V. (x) = 1-V (x) x <x<X 
ny nol a ny 
] x>xX 
ny 
The U;(Y) are dual. Then, if we define 
W; jy) = V,G)U5(Y), i= Ty...oMys j= Taessomy . 


it is easily observed that the function W, joy) has support c1(R; 5) and 
that the functions form a partition of unity for the planc, i.c., 
Hie hice = 1 forall (x,y) . 

These properties allow the construction of the interpolation function 
(1) to proceed easily since any point (x,y) is in at most four Ray » and 
the denominator of (1) is always = 1 , allowing us to write 


F(x,y) = Ps eet tags . 


We again emphasize that at most four terms in the sum are non-zero. 


«S. 
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The appropriate choice of x, and yj as well as ny and ny depend 
on the data as well as the choice of local interpolating functions Q; (xy) ; 
For this reason we defer discussion of the selection of these grid lines until 


after we discuss the choice of Q, 5 %-y) : 


2.2. Local Interpolation Functions 


The only restriction on the local interpolation functions Q, jy) are 
that they interpolate all data points in Rig and that they are defined for 
all (x,y) in Raj . Polynomials sometimes fail to satisfy these conditions. 
The use of optimal approximations in Sard corner spaces has been investigated 
[5], and for small numbers of data points, the approximations can be computed 
in straightforward fashion. One possible defect in such approximations is 
their lack of polynomial precision: even constants are not approximated exactly. 
With only a slight complication this can be overcome, since by a theorem of 
Barnhill and Gregory [3], the boolean sum operator B@L has the interpolation 
properties of B and the function precision of L . Here we are thinking of 
Bf as the optimal approximation in Bro .21 while Lf is the least squares 
fit by a linear function. 

The implemented version of the program embodies three options: (1) Use 
optimal approximations in Bro 2] as the local interpolation functions; (2) 
Use the least squares linear approximation instead of an interpolation function; 
and (3) Use the optimal approximation in Bro,2] boolean sum the least squares 
linear approximation. The second option yields a surface which in general does 
not interpolate the given data. The third option is achieved computationally 
as (BOL)f = (B +L - BL)f = Lf + B(f - Lf) . 

The use of the boolean sum has a desirable effect in that it removes much 


of the effect of linear transformations of the data on the overall approximation. 
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However, for complete consistency with respect to translation and change 


of the measure of distance, each rectangle 
X51 2X p47 ALY 5_ 7 Yg47] 
is transformed to the unit square for the optimal approximation. For these 


purposes, we take Xo = min Xp and x = Max Xp 5 and the dual in y. 


n+] 
The base point (a,b) is taken to be (0,0) for all i,j => 2 , while it is 
taken to be (1,1) for i=j=1, (0,1) for i22, j=1, and (1,0) 

for i=1, j 22. This yields lines of discontinuity in the second deriv- 
atives which are nowhere interior to the support regions for the weight func- 
tions Wy » thus assuring continuous second derivatives in the overall approx- 
imation. 

The overall approximation is invariant with respect to linear transfor- 
mations which leave the directions of the axes unchanged. Since lines of 
discontinuities in the third derivatives occur along horizontal and vertical 
lines the approximation is not invariant with respect to rotations. 

The points associated with Rij include all the points in the closure 
of Rig . Because approximation by a linear function requires at least three 
points, a parameter MINPTS , is used to assure that at least MINPTS points 
are selected for each Rij . If extra points are required, they are taken as 
the closest points in the sup norm, distance being measured after 

D5 ja XLV; Y 547] 
has been transformed onto [0,1]° . Presently MINPTS is set to three and 
this has been satisfactory. It is easily changed, if desired or necessary. For 
example, if some Rij has only three colinear points associated with it, the 


scheme will fail under options (2) or (3). Then one must either increase the 


value of MINPTS or use option (1). 


2.3. Selection of Grid Lines 

It is desirable to have automatic selection of grid lines, that is, 
values of x and iF . This should be accomplished in some manner which 
results in rectangles Rij which contain approximately equal numbers of 
points. For data which is poorly distributed this may not be possible. 
However, for somewhat uniformly distributed points the process we describe 
here works quite well. 

The selection of the grid lines is determined by one parameter, called 
NPPR, for “number of points per rectangle." The grid lines are then chosen 


so that there will be approximately NPPR points in each rectangle, R If 


i * 
there are additional points added to certain rectangles to make up MINPTS 
points the average may be higher. The average is, of course, dependent on 
the data set. 
Equal numbers of grid lines are chosen in each direction, that is le 
Because we want NPPR points per rectangle, each subrectangle 
(5 oa KY 5 9V jay) 


should have NPPR points. Thus we want to choose n, =n, so that 


X y 
(n, + 1)? . 7 NPPR =n _, the total number of data points. Thus, we take ny 


pl 


to be the nearest integer to (4n/nppR) '/2 eg 

Grid lines, that is x, and VF values, are now determined by choosing 
these values so that approximately n/(n, + 1) ee occur in each 
(Ks sXe 47] and each (V5 ¥ja1] . Specifically, let x, denote the values of 
Xy given in nondecreasing order, then X * x, » Where k_ is the integer 
nearest in/(n, +1) for i= 1,2,...9n, . The selection in y is dual. 
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3.0. Implementation 


The scheme is implemented in a set of subprograms, only one of which is 
normally referenced by the user. The hierarchy of subprograms is given in 


figure 1. A brief description of them, according to level, follows. 


LOCLIP CLOB 22 EVLB 22 


LEQT1P* 


Figure 1 


3.1.1. User's Program 

The user's program must provide the data points, (X,Y, oF)» k=1, 
.+.9M 5 aS well as the points x0; and y0; for the grid of points at which 
the interpolation function is to be evaluated. In addition the user's program 
must provide workspace arrays, IWK and WK, and an array FO for the returned 
function values. The amount of storage required in the arrays IWK and WK are 
not known a priori, but are estimated as follows. 

For IWK, approximately 4n(1 + 1/NPPR) locations are required. This 


has generally proven to be an overestimate. For WK, the storage required 


depends on the type of local approximation used and is approximately: 


*This is a Cholesky decomposition equation solver in the IMSL Library, which 
itself references several other subroutines from the IMSL Library. 
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for MODE = 1, 4(n + Yn/NPPR); 


for MODE = 2, 4(3n/NPPR + Yn/NPP®*; and 
for MODE = 3, 4(n + 3n/NPPR + Yn/NPPR) . 


As with the estimates for IWK, these estimates have proven to generally be 
overestimates. The precise number of storage locations required in each 
array are returned to the user's program by the principal subroutine. 

Under the usual option, the user specifies NPPR > 0 . If the user 
wishes to specify the grid lines, NPPR is set to zero, and additional input 
in the arrays IWK and WK is necessary. In the array IWK, the user specifies 
ny in IWK(1) and ny in IWK(2). The grid lines are specified in the array 
WK, according to the following. 

WK(1) is min Xp 

WK(2),...5WK(n, + 1) are the vertical grid values Xp oe eax in 
increasing order, 


WK(n, +2) is * Xp 


WK(n, + 3) ts wy Ye? 


WK(n, + 4) 4... 5WK(n, + ny + 3) are the horizontal grid values 

ce in increasing order, and 

WK(n, + My +4) is “ns Ye 
3.2.1. Subroutine LOB 22 

This subroutine provides the interface between the user and the set of 
subroutines which implement the method. Generally LOB 22 sets up storage 
locations in the arrays IWK and WK, determines parameters required by other 
subroutines, and calls subroutines to (1) generate the grid (if necessary), 
(2) determine the interpolation points for each rectangle, (3) compute coeffic- 
ients for the local interpolating functions, and finally, (4) to evaluate the 
interpolation function on the desired grid of points. 
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3.3.1. Subroutine GRID 

This subroutine selects the values of x and Y; in accordance with 
the discussion in Section 2.3. 
3.3.2. Subroutine LOCLIP 

This subroutine determines the local interpolation points for each Rij 
in accordance with the last paragraph of Section 2.2. 
3.3.3. Subroutine CLOB 22 

This subroutine computes the coefficients in the least squares plane 
(MODE = 2 or 3) and the coefficients in the optimal approximation (MODE = 1 
or 3) or each of the rectangles Ri . In the present implementation the 
IMSL Cholesky decomposition equation solver LEQTIP is used. This could be 
replaced at a facility where IMSL is not available, although according to 
IMSL policy, LEPTIP (and associated subroutines) can be used as part of this 
package at any facility. Because of the short single precision word length 
on the IBM 360/370 series computers, on which this program was implemented, 
the coefficients for the system of equations for the optimal approximation 
are generated in double precision. On computers with a longer word length, 
the double precision variables in this routine can be safely removed. Other 
double precision statements occur in EVLB 22 and function K, which must also 
be removed. 
3.3.4. Subroutine EVLB 22 

This subroutine evaluates the approximation (2) on the set of points 
(x0, -y0,) as specified by the user, and returns the values in FO. As noted 
above, the double precision variables in this subroutine should be removed 


on computers with longer word length than the IBM 360/370 series. 


mie 


3.4.1. Function K 
This function evaluates the representers for point evaluation functionals 
in Bro .2] . For evaluation at (u,v) , base point at (a,b) , the repre- 
senter, as a function of (s,t) [6] is 
K(a,b3u,v,S,t) = g,(asu.s)go(bsv,t) » where 
g,(asu,s) = Go(asu,s) = (s-u){9) + 1 + (u-a)(s-a) + (u-a) (s-a) (2) + (s-a) (3) 


for as<u_, and 
Jo(asu,s) = G,(-a3-u,-s) for u<a. 

The arguments of this function are all single precision, but because of 
the short word length of the ‘BM 360/370 computers , all calculations are per- 
formed in double precision, and the returned value is double precision. On 
computers with longer word lengths these calculations can be done in single 
precision. 

4.0. Examples 

The method has been applied to a number of sets of data with good results. 
Figures 3 - 5 show test surfaces and results of applying the method for each of 
the options for local approximations. The three surfaces are described by 


(C) F(x,y) = tanh (y - x) +1, 


2 
(S) F(xsy) = 3/2[cos(3/5(y - 1)) + 5/4]/(1 + (4) 1, and 
2 2 2 
(E) F(xyy) = 9(3/4 exp(=Ue3) 5 3) ) + exp(-) - ()) 


- & exp(-(x - 5)° = (y - 8)*) 
2 2 
+ F exp(-8) 7 tvs), - 
respectively. The 100 interpolation points were chosen at random within a unit 
square centered at (i,j) for i,j = 1,2,...,10 . The points are shown in 


Figure 2 as +'s, with the convex hull shown by dashed lines, while the square 
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(1,10)? » on which the resulting interpolation functions were evaluated, is 
given by the solid lines. The diagonal line shows the direction toward the 
viewing point. 

There does not appear to be a great deal of difference between the optimal 
approximation and the optimal approximation boolean sum least squares plane. 
Generally the latter option has slightly smaller errors and slightly less notice- 
able defects. Gross defects in the approximations can generally be traced to a 
lack of data in that particular part of the region. 

The effect of varying the parameter NPPR is shown in Figures 6 - 14. Some 
general observations are possible from this setof views. Most apparent is the 
fact that option (2), the least squares plane fit as the local approximation 
does not appear to lead to very good results. In general, however, the smaller 
value of NPPR gives better results, visually, and usually better accuracy, too. 

The choice of NPPR = 6 for options (1) and (3) appears to be a reasonable 
One. For surfaces with sharp gradients, as in Figure 6, it appears that local- 
izing the behavior as much as possible with a smaller value of NPPR is the best 
strategy. For smooth surfaces, such as in Figure 9 and 12 it appears the opposite 
is true, where NPPR = 8 seems to lead to the best results. 

The storage and timing results are given in Table ]. The storage refers 
to requirements of the two workspace arrays provided by the user. The timing 
is for calculation of the 1089 points generated for the plots. The program was 
run under the Fortran H compiler on the IBM 360 model 67 at the Naval Postgraduate 
School. Computation times are dependent on external factors and may vary from run 


to run. 
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Figure 2 
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NPPR MODE NIWK NWK Approx. Time (sec. 
4 1 437 375 7.9 
4 2 437 265 1.3 
4 3 437 618 8.2 
6 1 380 346 10.9 
6 2 380 165 1.1 
6 3 380 493 11.1 
8 1 351 328 12.9 
8 2 351 124 1.0 
8 3 351 436 13.1 
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Figure 5 (NPPR = 6) 
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Figure 11 (Mode = 3) 
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Appendix: Program Listings and Sample Output 
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